Output-input coupling in thermally fluctuating biomolecular machines 



Michal Kurzynski/'Q Mieczyslaw Torchala/'^ and Przemyslaw Chelminiak^ 

^Faculty of Physics, A. Mickiewicz University, Umultowska 85, 61-614 Poznan, Poland 
^ BioInfoBank Institute, Limanowskiego 24- A, 60-744 Poznan, Poland 
(Dated: June 27, 2011) 

Biological molecular machines are proteins that operate under isothermal conditions hence are re- 
ferred to as free energy transducers. They can be formally considered as enzymes that simultaneously 
catalyze two chemical reactions: the free energy-donating reaction and the free energy-accepting one. 
Most if not all biologically active proteins display a slow stochastic dynamics of transitions between 
a variety of conformational substates composing their native state. In the steady state, this dy- 
namics is characterized by mean first-passage times between transition substates of the catalyzed 
reactions. On taking advantage of the assumption that each reaction proceeds through a single pair 
(the gate) of conformational transition substates of the enzyme-substrates complex, analytical for- 
mulas were derived for the flux-force dependence of the both reactions, the respective stalling forces 
and the degree of coupling between the free energy-accepting (output) reaction flux and the free 
energy-donating (input) one. The theory is confronted with the results of random walk simulations 
on the 5-dimensional hypercube. The formal proof is given that in the case of reactions proceeding 
through single gates, the degree of coupling cannot exceed unity. As some experiments suggest such 
exceeding, looking for conditions of increasing the degree of coupling over unity challenges theory. 
Though no analytical formulas for models involving more transition substates are available, study 
simulations of random walks on several model networks indicate that the case of the degree of cou- 
pling value higher than one occurs in a natural way for scale-free tree-like networks. This supports a 
hypothesis that the protein conformational transition networks, like higher level biological networks: 
the proteome and the metabolome, have evolved in a process of self-organized criticality. 

PACS numbers: 05.70.Ln, 87.15. Ad, 87.15. Ak, 87.15.Hp, 87.15.Ya, 89.75,Hc 



I. INTRODUCTION 

An almost common conviction that biochemical pro- 
cesses can be interpreted in terms of the conventional 
chemical kinetics is based on an assumption that inter- 
nal dynamics of biomolecules is fast enough to ensure 
reaching the partial equilibrium state before each kinetic 
step [ij. However, at least a decade ago it was already 
clear that this assumption cannot be true, as besides fast 
vibrations the dynamics of biomolecules comprises also 
slower stochastic transitions between a variety of con- 
formational substates [2j-t6|. Research of biomolecular 
dynamics is being developed faster and faster. Today, 
even in the case of small, relatively rigid water-soluble 
proteins, one speaks about the 'native state ensemble' 
[7h14| rather than a single native state earlier identified 
with the protein tertiary structure, and for very small 
proteins or protein fragments trials to reconstruct the ac- 
tual networks of conformational transitions are realized 
[T5| - |2]| . As many as 30% native proteins is considered 
to be either completely or partly intrinsically disordered 
[23 - l23 | . Folding of such proteins takes place during the 
binding to targets 1251 [26| what is essential for molecu- 
lar recognition [2?!, 1281]. Koshland's concept of induced 
fit has been replaced by the 'conformational selection' 
concept [29l - [33 [. Allosteric regulation app ears to have 
dynamic rather than structural nature |29l437l |. 
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Because of the slow character of the conformational 
dynamics, both chemical and conformational transitions 
in an enzymatic protein have to be treated on an equal 
footing 6] and jointly described by a system of coupled 
master equations 

Pi{t) = ^[wii'Pi'it) - wi,ipi{t)] , (1) 
I' 

determining time variation of the occupation probabili- 
ties pi{t) of the individual protein's substates (Fig. [T]). 
In Eq. ([T]), Win is the transition probability per unit time 
from the substate / to I' and the dot denotes the time 
derivative. The conformational transition probabilities 
satisfy the detailed balance condition which, however, 
can be broken for the chemical transition probabilities 
controlled by concentrations of the enzyme substrates. 

In the closed reactor, a possibility of a subsequent 
chemical transformation to proceed before the confor- 
mational equilibrium have been reached in the actual 
chemical state results in the presence of a transient non- 
exponential stage of the process and in an essential dy- 
namical correction to the reaction rate constant describ- 
ing the following exponential stage [l|, [2, lol-fl^ . In the 
open reactor under stationary conditions (the concen- 
trations of reactants and products of the reaction kept 
constant), the general situation is more complex but for 
reactions gated by single conformational transition sub- 
states (Fig.[IJc)) a simple analytical theory was proposed 

[s^ . A consequence of the slow conformational tran- 
sition dynamics is that the steady-state kinetics, like the 
transient stage kinetics, cannot be described in terms of 
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FIG. 1. Development of kinetic schemes of a single enzymatic reaction R -o P. (a) Two-step Michaelis-Menten kinetics 
involving one enzyme-substrate intermediate M. (b) Three-step Haldene's kinetics involving two intermediates. Here, M' = 
ER and M" — EP. (c) Kinetics studied in Ref. where transitions between intermediates within E and M were expanded 
to quasi-continuous networks of conformational transitions described by parts of Eqs. ([1} and represented here by the gray 
boxes. The reactant and product binding-releasing reactions are assumed to be gated, i.e., they take place only in certain 
conformational substates represented here as black dots, (d) Simplified kinetic scheme considered in the present paper. Two 
reactant and product binding-releasing reactions and a kinetics of transitions within E are replaced by a single bimolecular 
reaction. All reactions are reversible; the arrows indicate directions assumed to be forward (the corresponding rate constants 
in text are written with the subscript -f). 



the usual rate constants. This possibility was suggested 
almost forty years ago by Blumenfeld [39]. Later on, we 
have shown that more adequate physical quantities that 
should be used are the mean first-pa ssag e times between 
distinguished transition substates d, [Sq . 

An application of the formalism to two coupled en- 
zymatic reactions was considered in the context of the 
free energy transduction in biological molecular machines 
[ssj . We understand the word 'machine' quite generally 
as denoting any physical system that enables two other 
physical systems to perform work one on another. Under 
isothermal conditions, performance of work is equivalent 
to a transduction of free energy. Thus, molecular ma- 
chines that operate under such conditions are referred to 
as free energy transducers [40j . 

From a theoretical point of view, it is convenient to 
treat all biomolecular machines, also pumps and motors, 
as chemo-chemical machines enzymes that simulta- 
neously catalyze two chemical reactions: the free energy- 
donating reaction and the free energy-accepting one. Un- 
der isothermal conditions, all chemical reactions proceed 
due to thermal fluctuations: a free energy needed for 
their realization is borrowed from the environment and 
then returned to it. In fact, the biological molecular ma- 
chines are Maxwell's demons: their mechanical or electri- 
cal elements are 'soft' and perform work at the expense 
of thermal fluctuations |4l| - t43| . Of course. Maxwell's de- 
mon can operate only out of equilibrium and it is a task 
of the free energy-donating reaction to secure such con- 
ditions. 

For the chemo-chemical machines, the degree of cou- 
pling, i.e., the ratio of the free energy-accepting (output) 
reaction flux to the free energy-donating (input) one was 
found to be also determined by the mean first-passage 
times between the conformational substates forming the 
reaction gates. As the mean first-passage times are not 
the quantities that could be directly determined in ex- 
periment, no experimental verification of the theory pre- 



sented in Ref. [s^ has been done as yet. The first goal of 
the present paper is to check the correctness of the theory 
by confronting it with results of Monte Carlo simulations 
performed on simple model networks of conformational 
transitions as well as to introduce quantities that could 
be determined experimentally. 

The essential motive of our studies is a trial to answer 
the intriguing question whether is it possible for the de- 
gree of coupling to have a value higher than unity. A 
dogma in the physical theory of, e.g., biological molec- 
ular motors is the assumption that for making a single 
step along its track the motor molecule has to hydrolyze 
at least one molecule of ATP Several years ago, this 
assumption has been questioned by a group of Japanese 
biophysicists from the Yanagida laboratory who, joining 
a specific nanometry technique with the microscopy fiu- 
orescence spectroscopy, shown that the myosin II head 
can make several steps along the actin filament per ATP 
molecule hydrolyzed |45l447| . This observation has been 
confirmed by some other laboratories [48?], also for the 
dynein [49ml| moving along the microtubules. In Refs. 
[3 a ) and basing on approximations carried too far, we 
suggested that the degree of coupling can exceed unity 
already for reactions proceeding through single pairs of 
transition substates. Here, we formally prove that it is 
not the case and show that the latter possibility real- 
izes a natural way for the scale-free tree-like models of 
conformational transition networks with the output gate 
extended to more conformational substates. 



II. GENERALIZED KINETIC SCHEME OF 
CHEMO-CHEMICAL MACHINE ACTION 

The principle of action of the chemo-chemical machine 
is simple (40| . It is a protein enzyme that catalyzes simul- 
taneously two chemical reactions (Fig. (^fa)). Separately, 
each reaction takes place in the direction determined by 
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FIG. 2. (a) Principle of chemochemical free energy transduction. Due to proceeding on the same enzyme, reaction Ri •<->■ Pi 
drives reaction R2 -f-*- P2 against its conjugate force determined by steady state concentrations of the reactant and the product, 
(b) Assumption of a possible short circuit or slippage of the input vs. output reaction, (c) Assumption of both the free 
energy-donating and the free energy-accepting reaction to participate in a kinetic scheme like the one shown in Fig. [TJd). 



the second law of thermodynamics, i.e., the condition 
that energy dissipated, determined by the product of 
flux and force, is positive. However, if both reactions 
take place simultaneously in a common cycle, they must 
proceed in the same direction and the direction of the 
first reaction can force a change of direction of the sec- 
ond. As a consequence, the first reaction transfers a part 
of its free energy recovered from dissipation performing 
work on the second reaction. 

In formal terms, the chemo-chemical machine couples 
two unimolecular reactions: the free energy-donating re- 
action Ri -f-j" Pi and the free energy- accepting reaction 
R2 -H- P2. Bimolecular reactions can be considered as 
effective unimolecular reactions on assuming a constant 
concentration of one of the reagents, e.g. ADP in the case 
of ATP hydrolysis. The input and output fluxes Ji (i = 
1 and 2, respectively) and^the conjugate thermodynamic 
forces Ai are defined as [4C| 



J,. = 



[E]o 



and 



PA, = In K, 



K,, = 



[P. 



oq 



eq 



(2) 



(3) 



Here, symbols of the chemical compounds in square 
brackets denote the molar concentrations in the steady 
state (no superscript) or in the equilibrium (the super- 
script eq). [E]o is the total concentration of the en- 
zyme and (3 is proportional to the reciprocal tempera- 
ture, /3 = (fceT)"^, where fee is the Boltzmann constant. 
The thermodynamic forces measure the distance from the 
equilibrium at which they vanish. The flux-force depen- 
dence is one-to-one only if some constraints are put on 
the concentrations [R^] and [Pi] for each i. There are 
two possibilities. Either the concentration of one species, 
say Ri, in the open reactor under consideration is kept 
constant: 

[R,] = const. = [Rif'i (4) 
or is such the total concentration of the enzyme substrate: 
[Ri] + [P,] = const. = [R],o . (5) 



The free energy transduction is realized if the product 
J2A2, representing the output power, is negative. The 
efficiency of the machine is the ratio 
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(6) 



of the output power to the input power. In general, the 
degree of coupling 



e = J2/J1 



(7) 



being itself a function of the forces Ai and A2, can be 
both positive and negative. 

Usually, the assumption of tight coupling between the 
both reactions is made (Fig. [H^a)). It states that the 
flux of the first reaction equals the flux of the second, 
Ji — J2 thus e = 1. However, an additional reaction 
can take place between the two states M' and M" of the 
enzyme-substrates complex (Fig. [1Kb)). The latter re- 
action can be considered either as a short circuit, the 
non-productive realization of the first reaction not driv- 
ing the second reaction, or a slippage, the realization 
of the second reaction in the direction dictated by its 
conjugate force. The short circuit and the slippage are 
more favorable thermodynamically, so for the free energy 
transduction still to take place, the remaining reactions 
should be more favorable from the kinetic standpoint. 

The multiconformational counterpart of the scheme in 
Fig. [21(b) is shown in Fig.[2Kc). Here, like in the scheme in 
Fig. [TJ^d) , a network of conformational transitions within 
the enzyme-substrates complex is represented by the gray 
box and the assumption of gating by single pairs of con- 
formational transition substates is made. The input and 
the output reaction fiuxes are determined by station- 
ary occupation probabilities of the gating transition sub- 
states. These can be calculated with the help of a tech- 
nique of summing up the directional diagrams proposed 
by Terell L. Hill i4Q] wh o formalized an old idea of Gustav 
Kirchhoff. In Ref. [38|, we shown how various classes of 
directional diagrams for the networks of conformational 
substates are related to the mean first-passage times be- 
tween distinguished substates. Hence, the fluxes are to 
be expressed in terms of appropriate mean flrst-passage 
times. In the next Section, we quote the most important 
results of Ref. [11] in the essentially changed notation 
facilitating their direct experimental verification. 
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III. SINGLY GATED REACTIONS: 
THEORETICAL RESULTS 

For all the schemes shown in Fig. [21 the flux-force de- 
pendence for the two coupled reactions has a general 
functional form [33 



1 



-fi(Ai--Af) 



J 



(8) 



The parameters J+,;, J_i, Joi and Af depend on the 
other force and are determined by a particular kinetic 
scheme. Af" have the meaning of stalling forces for which 
the fluxes Ji vanish: Ji{Af) = 0. The dependence Ji{Ai) 
is strictly increasing with an inflection point, determined 
by Joi, and two asymptotes, J+i and J-i (compare Fig.[4| 
further on). The asymptotic fluxes J+i and J_i display 
the Michaelis-Menten dependence on the substrate con- 
centrations: 



K 



k+i [Ri] 



(9) 



Because of high complexity, we refrain from giving any 
formulas for the turnover numbers k±i and the apparent 
dissociation constants K±i. 

Without the loss of generality we assume that both Ji 
and Ai are positive. However, J2 can be either positive 
or negative. Accordingly, the stalling force should 
be negative or positive. Within the range of free en- 
ergy transduction, A|* < yl2 < or < yl2 < Af, the 



flux-force dependence 72(^2) can be convex, concave or 
involving an inflection point as well. The linear Onsager 
approximation is in general not applicable. 

The degree of coupling dependence on the forces Ai 
and A2 has a functional form 



and the stalling force 



/3Af =ln- 



-PA 



^+Wi{Ai) 
l + Wi{Ai) 



(10) 



(11) 



The expression for Af is to be obtained after replacing 
the index 1 with the index 2 and vice versa. The quan- 
tities Wi{Ai) are measures of the slippage. The lower 
slippage the lower value of the corresponding Wi. For 
no slippage we simply have A|* = — Ai, Af = — A2, and 
e = 1 (the both fluxes Ji and J2 are driven by the same 
resultant force Ai -I- A2). For finite slippage, Af reach 
the maximum, negative or positive values in the asymp- 
totic limits of sufficiently high values of the other forces 
(compare Fig. [Bj further on). 

For the simple scheme shown in Fig. [IJb), a direct 
calculation results in 

W^i(Ai) = A:o+ri(Ai), ^^2(^2) = ^0^x2(^2) , (12) 

whereas for the scheme in Fig. [2l^c) , the summation over 
diagrams gives [ssj 



J 



Wi{Ai) 



tm{1"^{1', 2'}) + rM(l^o{l^ 2"})e-P^- + tMi) 
TM(l'o{l",2'})-rM(l'o{l",2"}) 



(13) 



and the adequate for ^^2(^-2) after replacing 1 with 2 
and vice versa (note the 180° rotational symmetry of the 
kinetic scheme in Fig. [He)). The quantities Ti{Ai) in 
Eqs. (fT2]) and (|T3| are mean transition times in the for- 
ward direction outside the enzyme-substrates complex M. 
For the outside transition of the form shown in Fig. [IJc), 
this time is given by the expression [sl] 



r{A) 



ik'-] 



(14) 



P°q(M)i^rE(0'oO") 



-I3A 



(te(0' f+ 0") is the sum of the forward and the reverse 
mean first-passage times between the distinguished sub- 
states in E) . For the outside transition of the form shown 
in Fig. [Hd), it is simplified to 



riA) = (fc+[R])- 



(15) 



(16) 



I 

constants defined as 

k'l^K'lpl?,{M)/P^^{M), 
k'_ = k'_p°?(M)/P'='1(M), 
fc+ = ^+pS?,(M)/P°^(M). 

pJ^'^(M) denotes the equilibrium occupation probability 
of the transition conformational substate I in M and 
P'='3(M) denotes the equilibrium occupation probability 
of the whole enzyme-substrate complex M. Notation of 
the transition probabilities between reaction transition 
substates per unit time k", k' and k is explained in Fig.[T] 
The quantities tm in Eq. (jl3p denote the sums 



Above, we introduced the transition state theory rate 



TM{lo^{lJ'})=TM{lo^{lJ'}) + rM{{l,l'}'^lo) (17) 

of the forward and reverse mean first-passage times that 
occur in the summation formula for the mean first- 
passage time from Iq to / in the network symbolized by 
the gray box M: 

Tuilo-^l) = rM{{lo,l'}^l) +tm{Io^{1,1'}) (18) 
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for arbitrary V . Eq. (fT8|) is a generalization of the obvious 
summation formula for one-dimensional diffusion: 

t{Io^I)=t{1^^1')+t{1'^1), (19) 

V lying between and L The quantity tm(^o {^^'}) 
has a direct meaning of the mean first-passage time from 
lo to I or I'. The interpretation of tm{{Io,1'} ^l) is more 
troublesome but it can be always treated as a completion 
of tm{Io^{1,1'}) to tm(^o— ^0- On doing it we find two 
alternative relations 

tm{Io^{IJ'})^tm{Io-^{1,1'}) (20) 

= TMil0'^{l'j}) 

-tm{1'-^{Io,1}) + tm{1'-^Io) ■ 



IV. SINGLY GATED REACTIONS: 
COMPARISON WITH MONTE CARLO 
SIMULATIONS 

The quantities tmHq-^ {I, I'}) occurring in Eq. (fT3|) for 
Wi(Ai) and the adequate for ^2(^2) can be considered 
as six independent parameters of the theory to be fitted 
in future experiments. However, the mean first-passage 
times occurring in Eqs. (pUj) are not the quantities that 
could be directly determined experimentally. The choice 
of an appropriate network that models the interior of the 
gray box in Fig. [2fc) and positions of the input and out- 
put gates is a question of the statement of a more or less 
reasonable hypothesis. The simplest was stated seventy 
years ago by Kramers who assumed that the slowly vary- 
ing intramolecular substates lie along a one-dimensional 
'reaction coordinate' 0, H^]- This way of the reasoning 
was continued in the theory of molecular motors where 
the reaction coordinate was identified with the position 
of the motor along its track (io Hil lssI, Isi] . 

More complex modeling can base on statistical anal- 
yses of time series found in molecular dynamics simula- 
tions [T5| - |2]| or single-molecule experiments [55l - [6^ . Un- 
fortunately, the present-day knowledge in this matter is 
still rather poor, so we decided first to test our theory by 
resorting to Monte Carlo simulations of random walks 
on simple, but not quite real networks. Various networks 
differ one from other by the geometry of links that deter- 
mines an entropic contribution to the kinetics, and the 
variety and asymmetry of the transition probabilities ww 
that, following the detailed balance condition, determine 
an energetic contribution. 

In the beginning, for more detailed studies we chose 
the most regular isoenergetic network, the n-dimensional 
hypercube, the vertexes of which are labeled by sequences 
of the bits (si, S2, . . . , Sn), Si = 0, 1, and all possible tran- 
sitions are related to the change of one bit and have the 
same probability w. The distance between vertexes is 
determined by the minimum number of edges a random 
walker has to pass in a walk between these vertexes. The 



such determined distance equals the number of the neces- 
sary bit changes. In the n-dimensional hypercube, there 
are N = 2" vertexes, each vertex has n neighbors and 
no boundary conditions are necessary. For a reasonable 
time of simulations we assumed the dimension n — 5. 

We chose the input and the output gates so as to make 
the free energy transduction the most effective. It takes 
place when moduli of both the degree of coupling ([25)1 
and the stalling force (fTTj) are maximum, i.e., the values 
of W2{0) and VFi(oo) are minimum. A detailed analysis 
indicated that it holds if the pairs of sites 1' and 2" as well 
as 2' and 1" lie at the closest, at the distance equal to 1. 
On the contrary, the diameters of the input and output 
gates should be larger. We chose them the largest, equal 
to 5 so as to lie along the diagonals of the hypercube. 

The such determined geometry of the gates is unique 
and, moreover, only the values of three types of the mean 
first-passage times have to be known, enabling one to 
calculate, following Eq. (PO)) . all the quantities occurring 
in Eq. ^ for Wi{Ai) and the adequate for 14^2(^2)- 
These are the mean first-passage times of the type 



tm(1"-^{1',2'}) = 16, 
tm(1' ^ {1", 2'}) = 80/3 w 26.667 , 
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FIG. 3. Simulated time course of the net number of the input 
(Ri -f-^ Pi) and the output (R2 •<->■ P2) external transitions for 
the 5-dimensional hypercube with gates and parameters de- 
scribed in text, (a) Snapshots made every step, (b) Snapshots 
made every 10^ steps. 
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FIG. 4. Flux-force dependence J2{A2) for the model and pa- 
rameters described in text. The black dots denote results 
of the Monte Carlo simulations and the continuous line rep- 
resents the fit to Eq. The free energy transduction is 
realized in the range —0.48 < l3A2 < 0. 

tm(1'^1") = 128/3 « 42.667 

(note that for the hypercube, the mean first-passage 
times depend only on the distances between the ini- 
tial, final and the intermediate, if any, states). All the 
quoted values, counted in the number of the random walk 
steps with the transition probability between the nearest 
neighbors w = l/n = 1/5, were determined by simple 
though tedious combinatorics and checked in numerical 
simulations. 

We assumed the simplest, constant reciprocal forward 
external transition times given by Eq. (|15p with the use 
of Eq. Q and chose rf ^ = 50 w/N and t^^ ^ 30 w/N, 
counted in the reciprocal time units w. For w — l/n — 
1/5 and the equilibrium occupation probability of the 
each lattice vertex 1/N — 1/2" — 1/32, the times 

Ti = 3.200, T2 = 5.333 

are one order of the magnitude shorter than the maxi- 
mum mean first-passage time tm(1' — 1") being a mea- 
sure of the intramolecular relaxation time. Thus, the 
both reactions 1 and 2 are controlled, though not com- 
pletely, by the intramolecular dynamics [l| . The recipro- 
cal reverse external transition times equal t^^ and 
multiplied by the detailed equilibrium condition-breaking 
exponents exp(— and exp(— /3A2), respectively. We 
chose /3Ai — 10 which is a physically reasonable condi- 
tion of the free energy donating reaction 1 to proceed 
sufficiently far from the equilibrium. 

The reciprocal external transition times multiplied by 
the equilibrium occupation probability l/N oi the tran- 
sition gates determine the external transition rates. In 
actual simulations, to preserve equal probabilities of the 
forward and the reverse internal transitions in the pres- 
ence of additional external transitions, w had to be cho- 
sen much lower than l/n and a high probability of wait- 
ing at each but one vertex had to be added. 
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FIG. 5. Dependence of the degree of coupling e on the output 
force A2 for the model and parameters described in text. The 
black dots denote results of the Monte Carlo simulations and 
the continuous line is calculated following Eqs. (|10p . (|13p and 
the adequate for 14^2(^42). No fit has been performed between 
experiment and the theory as the latter has no free parame- 
ters. The present figure is similar to Fig.|4]as for the very high 
value of the input force f5Ai = 10 the input flux Ji remains 
almost constant. 

A typical result of a simulation of the time course of 
the net number of external transitions through the input 
and the output gates is shown in Fig. [3lja). It is clearly 
seen that even for such the small lattice studied, consist- 
ing of 32 vertexes, large fluctuations make determination 
of the input and the output fluxes in 10^ iteration steps 
impossible. Only the increase of the number of the iter- 
ation steps to 10^ (figure [3ljb)) enables one to determine 
the fluxes with the error lower than 0.3%. Dots in Fig. |4] 
present the such determined values of the output flux J2 
as a function of the conjugate force A2. This figure also 
shows that the results of our Monte Carlo experiment fit 
very weU the theoretical prediction, Eq. ([8]). 

The theory presented in the previous section, in par- 
ticular its main Eqs. (fTO| . (fTTI) and ([13]), does not use any 
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FIG. 6. Dependence of the stalling force on the input 
force Ai for the model and parameters described in text. The 
black dots denote results of the Monte Carlo simulations and 
the continuous line is calculated following Eqs. ([TT|l and ([T3)l. 
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approximation and is exact. For simple networks of con- 
formational transitions like the one considered here, we 
know values of appropriate mean first-passage times and 
can compute directly the slippage functions Wi{Ai), thus 
the force dependences of the degree of coupling pO)) and 
the stalling force (fTTj). Figs. [5] and [6] show confronta- 
tions of the such obtained dependences with results of 
the Monte Carlo simulations that in the present context 
can be treated as experimental data. On taking into ac- 
count that in Figs. [S] and IHl no fitting procedures were 
applied, the agreement is excellent, but it only points to 
the correctness of conditions at which our Monte Carlo 
simulations were performed. Being sure of such correct- 
ness, we can use similar simulations for testing the qual- 
ity of various approximations we have to apply for actual 
networks. In particular, in two last sections we study by 
such methods effects of extending the gates to include 
many conformational transition substates. 



V. FOR SINGLY GATED REACTIONS THE 
DEGREE OF COUPLING IS LOWER THAN 
UNITY 

It is not simply to evaluate the range of values that 
the degree of coupling (|10l) assumes within the range of 
I3A2 corresponding to the free energy transduction. In 
Refs. [Ill and [l| , basing on inaccurate evaluations of the 
denominators in the expressions for Wi(v4i) and ^2(^2), 
we suggested that in the case of singly gated reactions, 
there are possibilities of the degree of coupling modulus 
having a value higher than unity. Here, we present a 
formal proof that this modulus cannot exceed unity. The 
proof makes use of two relations. 

From the symmetry basing the two equalities (P(7| , the 
first relation follows: 

tm{Io^{1,1'}) -ruil^^ilj"}) (21) 
^tm{1^{Iq,1"})-tm{1^{IoJ'})- 

It secures the expression under the logarithm in Eq. ((TT]) 
and the adequate for (3Af to be positive irrespectively 
of the sign of Wi . The second relation, not equivalent to 
(ICT , is to be derived from an equation analogous to 
involving two intermediate nodes: 

tm{Io^{1,1'}) -tm{Io^{1,1"}) (22) 

= TM{l'^{l",lo})~TMil'^{l"l}). 

From this relation, it follows that the denominators in 
M^i(Ai) and ^2(^2) equal each other. A consequence 
is that both Wi{Ai) are always of the same sign, either 
positive or negative. 

The highest value of the degree of coupling modulus in 
the free energy transduction region is for (3A2 — 0. Then, 
Eq. PU)) is simplified to 



independently of /3Ai, and 

/3^f(0)=0. (24) 

If the both Wi are positive than the stalling force f3Af 
given by Eq. is negative and the free energy trans- 
duction takes place in the region pAf < (3A2 < 0. The 
highest value of the degree of coupling ([23l) is positive 
but always lower than unity. 

One could expect higher value of the degree of cou- 
pling modulus if the both Wi were negative, i.e., if the 
denominator in (jl3|) was negative, which corresponded to 
interchange 2' with 2" thus J2 with — J2. However, it is 
not the case. The condition of the expression ((23|) to be 
lower than minus unity is 

W2{0r^ <-l/2. (25) 

On substituting the explicit expression for ^2(0) and 
taking advantage of the relation (|21|) this inequality can 
be rewritten as 

tm(2'o {1', 2"}) + rM(2"o {1", 2'}) + T2(0) < . (26) 

Of course, neither tm nor T2 can be negative, so the 
modulus of the coupling ratio can never be higher than 
unity. In consequence, if the both Wi are negative, the 
stalling force pAf given by (fTTj) is positive and the free 
energy transduction also takes place, now in the region 
< I3A2 < (3Af, however with the negative coupling 
ratio of the modulus also lower than 1 . 



VI. THE DEGREE OF COUPLING CAN BE 
HIGHER THAN UNITY. CASE OF SEVERAL 
SUCCEEDING OUTPUT GATES 

We proved the theorem that the value of the degree of 
coupling should be lower or at the most equal to unity, 
but only in the case when the input and output reactions 
proceed through single pairs of conformational transition 
substates. It is reasonable to suppose that a possibility of 
higher degree of coupling is realized if the output gate is 
extended to two or more pairs of the transition substates. 
Indeed, in Fig. [T^a) a scheme is shown with one input 
gate (l"a, I'a) and two succeeding output gates (2"a, 2'a) 
and {2"b, 2'b) closed in the common cycle. It is obvious 
that the degree of couphng for such scheme is e = 2. 
Similar reasoning has been proposed in order to explain 
multiple stepping of the myosin molecule along the actin 
filament [46|. One can imagine an incorporation of a 
system of additional transitions and an increase of the 
number of the transition substates what was for the first 
time considered by Terada and coworkers [63] . 

Unfortunately, even in the case of only two output 
gates the analytical formulas are so complex and not 
transparent that serious approximations are needed to 
be made from the very beginning. Being not able to for- 
mulate presently such approximations, we decided to ap- 
ply computer experiment for a preliminary study of the 
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problem. We performed Monte Carlo simulations start- 
ing from the 5-dimensional hypercube. For the most op- 
timal geometry of one input and two output gates we 
obtained the degree of coupling e(0) not higher than 
0.362. Also simulations on the 9-dimensional hypercube 
with 512 states were not successful. We suppose that 
the steady state fluxes on the isoenergetic hypercube of 
arbitrary dimension and with arbitrary number of gates 
are always equivalent to steady state fluxes on some non- 
isoenergetic (weighted) network with various probabili- 
ties of substate occupations and transitions with the sin- 
gle input and output gates. Such networks are described 
by the theory given above and all the discussion already 
performed applies to them. 

Still restricting ourselves to isoenergetic networks, we 
considered systems with more complex topology. We 
pointed our attention toward the networks with bot- 
tlenecks or dead ends, the diffusion on which displays 
long-time tiles jll]. The networks that model the ac- 
tual conformational dynamics ofproteins should display 
a hierarchy of relaxation times |2!-l5| . As an example of 
a network with a hierarchy of bottlenecks we considered 
the Sierpinski gasket (Fig. IH^a)) and as an example of a 
network with a hierarchy of dead ends, the Bethe lattice 
(Fig.mjb)). We assumed the values of the external transi- 
tion parameters as in the former section: = 50 w/N, 
T^^ = iOw/N and /3Ai = 10. For the Sierpinski gasket 
of the fourth order with boundary conditions and the sys- 
tem of gates shown in Fig. HJa), we got e(0) = 1.27. For 
the Bethe lattice with the five shells and the system of 
gates shown in Fig.|51[b), we got e(0) = 1.19. We conclude 
that for protein machines with the stochastic dynamics 
described by an appropriate network of conformational 
transitions, the degree of coupling can in principle be 
higher than unity. 

The geometry of gates shown in Fig. [5] was chosen 
with a bias against unfavorable short circuits or slippages 
and, simultaneously, long wandering between transitions 
through the successive gates. The goal was achieved in 
the evidently artificial way due to entropic obstacles and 
shortcuts. However, no obstacles and shortcuts could 



(a) 



(b) 




R 



FIG. 7. (a) Extension of the kinetic scheme in Fig. [JJ^a) to two 
output gates. Obligatory transitions are drawn by arrows. If 
no other transitions are realized, the degree of coupling equals 
two. Otherwise, it is lower than two but possibly higher than 
one. (b) Generalization of the scheme in Fig. [2{c) to a quasi- 
continuum of gates. The rate of the external transition within 
each gate can be different so that such kinds of models are 
referred to as the ones with 'fluctuating barriers' [Q]- 




FIG. 8. (a) Sierpinski gasket of the fourth order. Each vertex 
has four neighbors. The assumed boundary conditions allow 
additional transitions between the three vertexes of the largest 
triangle, (b) Bethe lattice with the flve shells. Each vertex 
besides the most external ones has three neighbors. For the 
most external vertexes, reflexive boundary conditions are as- 
sumed. In both networks, the geometry of the distinguished 
transition substates is shown. 



be needed if non-isoenergetic (weighted) networks were 
considered, with variable and appropriately chosen tran- 
sition probabilities. In this way models with fiuctuating 
barriers are to be obtained Q , symbolically presented in 
Fig. [Tl^b). Looking for simple model networks with the 
controllable and higher than unity degree of coupling that 
elucidates and, possibly, predicts the action of 'biomolec- 
ular gears' is certainly an important task both for the- 
oreticians and experimentalists. In the following, last 
Section, we present one possible proposal. 



VII. MORE NATURAL MECHANISM OF 
INCREASE THE DEGREE OF COUPLING IS 
OFFERED BY SCALE-FREE TREE-LIKE 
MODELS 

Since the formulation by Bak and Sneppen a cellu- 
lar automaton model of the Eldredge and Gould punc- 
tuated equilibrium j65j], the biological evolution is more 
and more often considered as a self-organized criticality 
phenomenon [6^, [67| . An evolving network model of self- 
organized criticality was proposed by Barabasi and Al- 
bert [6^ |69[ . Soon it appeared that two networks of the 
systems biology: the proteome and the metabolome have. 
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to a good approximation, not only the scale-free struc- 
ture like the Barabasi- Albert networks [t^, [zll but dis- 
play also a hierarchically modular, i.e., self-similar (frac- 
tal) organization [tI, [TJ]. An evolutionary mechanism 
has been proposed to elucidate the scale-free character 
of the proteome thus the metabolome [zl] ■ 

There are premises that also the conformational tran- 
sition networks in proteins are both scale-free and hi- 
erarchically modular. The former feature is suggested 
by results of molecular dynamics simulations for small 
atomic clusters [tI, 76 | and by a specific spatial orga- 
nization of proteins 77 , [tI] . The latter has been shown 



already in the pioneer papers from the Hans Frauenfelder 
laboratory Q and confirmed in early molecular dynamics 
simulations for the very proteins [3|-l5|- Thus, a hypoth- 
esis sounds reasonably, that also the protein conforma- 
tional transition networks have evolved in the process of 
self-organized criticality. 

However, evolutionary speculations above are quali- 
fied. The evolving scale-free Barabasi-Albert networks 
have not fractal but rather small- world character [69| . 
And, indeed, such a character was also suggested for 
both the proteome [79| and the conformational transition 
network [75l - l78l |. Only recently, an apparent contradic- 
tion between fractality and small- worldness have been ex- 
plained by application of the renormalization group tech- 
nique |80| . It appears that a network can be fractal in a 
small length-scale, simultaneously having the small- world 
features in the large length-scale and this is the case of 
the proteome and, probably, the protein conformational 
transition networks. 

The topological structure of the flow (of probability, 
metabolites, energy or information) through a network is 
characterized by a spatial spanning tree composed of the 
most conducting links not involved in cycles. It is referred 
to as the skeleton [8l[ or the backbone [82j of the network, 
all the rejected links being considered as shortcuts. The 




FIG. 9. Example of the Barabasi-Albert tree with 100 nodes 
and the dynamics described in text. The equilibrium occupa- 
tion probabilities of the nodes are distinguished by the size of 
the dot and its distance from the center of the circle formed 
out of the highest free energy nodes with one link. The input 
and the output reaction transition substates are shown. Note 
that there are seven output transition substates 2". 



skeleton of the scale-free network is also scale-free but 
the skeleton of the self-similar network needs not be self- 
similar. Here a criticality feature appears important that 
denotes the presence of a plateau equal to unity in the 
mean branching number dependence on the distance from 
the skeleton root. The critical skeletons can be completed 
to self-similar scale-free networks and such is the case of 
the proteome [Sll, Is^l • 

We state the hypothesis that such is also the case of 
the protein conformational transition network and that 
the plateau in the mean branching number versus the dis- 
tance dependence corresponds to the length-scale range 
where the original network displays the fractal properties. 
Diffusion on the fractal networks is characterized by the 
power-low first passage time distribution and a range of 
such distribution was found in Monte Carlo simulations 
of the random walk on large Barabasi-Albert trees [84|. 
In the corresponding length-scale range the mean branch- 
ing number tends to unity, i.e., the tree is in a sense 
equivalent to the one-dimensional chain transmitting the 
probability flow the fastest the possible. The latter prop- 
erty can be considered optimal for biological networks 
what justifies choosing the Barabasi-Albert trees as the 
object of our further considerations. 

Assuming the preferential attachment rule we can con- 
struct the Barabasi-Albert trees starting from a single 
node and adding one node in each construction step |69| . 
In Fig. 9, an example of such a tree is shown, obtained 
after 99 construction steps, thus having N — 100 nodes. 
To provide the network with a stochastic dynamics de- 
scribed by Eq. ([TJ , we assume the probability of changing 
a node to any of its neighbors to be the same in each ran- 
dom walk step. Consequently, the transition probability 
from the node I to the neighboring node /' 

Win = , 

where fc; is the number of links (the degree) of the node 
I. The network with such a dynamics cannot be isoen- 
ergetic and following the detailed balance principle the 
equilibrium occupation probability of the node / 



In the scheme presented in Fig. 9, the equilibrium occu- 
pation probabilities, i.e., the values of the corresponding 
free energies, are distinguished by the size of the dot and 
its distance from the center of the circle formed out of 
the highest free energy nodes with one link. 

As in Sections IV and VI we choose the simplest, con- 
stant mean forward external transition times given by 
Eq. with the use of Eq. ^ and assume their values 
to be much shorter than the longest mean first-passage 
time of the internal transitions. We choose ti — T2 — 2, 
which makes the external reactions almost completely 
controlled by the dynamics of the enzyme-substrate com- 
plex. The exit probabilities from the input {i — \) 
and the output (z = 2) reaction transition substates in 
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the forward direction (doubly primed) equal the recip- 
rocal product of the external transition time and the 
equilibrium occupation probability of that substates, cf. 
Eqs. (HH): 

The exit probabilities from the transition substates in the 
reverse direction (singly primed) are similar, but must 
be multiplied by the factor breaking the detailed balance 
symmetry, determined by the external forces: 

As in Sections IV and VI we choose f3Ai — 10 which 
makes the exit probability from the transition substate 
1' negligible. The exit probabilities from all the remain- 
ing transition states are much higher than the internal 
transition probabilities to the neighboring substates. For 
securing the sum of all transition probabilities from a 
given node to be one in the actual simulation step, all 
the discussed transition probabilities are appropriately 
renormalized and a high probability of waiting at each 
but one node is added. 

From what was told above it follows that both the in- 
ternal transition probabilities and the exit probability 
from a given transition substate are inversely propor- 
tional to its degree fc;. This fact implies a strategy of 
choosing the output gates. For the resultant output flux 
J2 at the zero external force {I3A2 — 0) to be the highest, 
the forward reaction transition substate 2" should have 
the lowest degree and be realized many times whereas 
the reverse reaction transition substate 2' should have a 
higher degee and be realized singly. As concerns the for- 
ward reaction transition substate 1", it should lie close 
to 2' and, because under the pumping conditions in the 
steady state its occupation decreases essentially, its equi- 
librium occupation should be the highest, i.e., it should 
be the main hub. The reverse reaction transition substate 
1' need not be highly occupied but it should lie close to 
any of the substates 2". The gates chosen in accordance 
with this strategy are shown in Fig. 9. 

In Figs. 10 (a) and (b), typical results of Monte Carlo 
simulations of the net number of external transitions 
through the input and the output gates, respectively, are 
shown. It is worth pointing attention to the 'devil's stair- 
case' form of the input flux in the range of medium tran- 
sition times, what is characteristic for diffusion on fractal 
networks. The absence of longer transition times results 
from the small-worldness effects discussed above as well 
as the boundary conditions. The absence of shorter tran- 
sition times results from the finite distance between the 
starting and the ending node as well as admitting addi- 
tional external transitions through the output gates. 

In the output flux, numerous reverse transitions are 
observed. Nevertheless, there are more forward transi- 
tions on the average per one input flux transition. In 
Fig. 10 (c), both fluxes are presented in much longer time 
scale from which the value of the output-input degree of 
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FIG. 10. Simulated time course of the net number of the input 
and the output external transitions for the model presented in 
Fig. [9] and parameters described in text (/3^2 = 0). (a) Snap- 
shots made every step for the input flux n\. (b) Snapshots 
made every step for the output flux 722. Many succeeding 
forward and reverse transitions are seen as broadened lines, 
(c) Snapshots for ni and 712 made every 10^ steps. 
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PA, 

FIG. 11. Dependence of the degree of coupling e on the output 
force A2 for the model presented in Fig. [9] and parameters 
described in text. The black dots denote results of the Monte 
Carlo simulations restricted to the free energy transduction 
region and the continuous line represents the fit to Eq. ((8|. 



coupling can be easy evaluated for e(0) = 11.0. The main 
reason for such a high value of e is a large representation 
of medium transition times in the input flux confronted 
with an approximately exponential distribution of tran- 
sition times in the output flux. 

Fig. 11 shows how the degree of coupling e decreases 
with the output force A2. For the model and the pa- 
rameters assumed, the stalling force can be evaluated for 
— 1.4 k^T units. It is worth noting an excellent fit to 
Eq. ([8|) what, on taking into account constancy of the 
input flux Ji within the region considered, suggests that 
the formula ^ remains applicable for any biomolecular 
machine. 



VIII. CONCLUDING REMARKS 

Only in a few recent years, some trials have been under- 
taken to determine conformational transition networks in 
native proteins. It is why in the present paper we re- 
stricted our attention to model networks. Our goal was 
to calculate and simulate the degree of coupling between 
the free energy-accepting and the free energy-donating 
reaction flux in protein molecular machines. Exact the- 
oretical formulas were possible to be obtained only for 
reactions proceeding through single pairs (the gates) of 
conformational transition substates. The theory predicts 
the value of the degree of coupling not exceeding unity. 
However, in Monte Carlo simulations on simple scale- free 
tree-like networks we shown that on increasing the num- 
ber of the output gates one can easily obtain the degree 
of coupling much higher than unity. In other words it 



means that 'biomolecular gears' are possible. 

Nevertheless, the degree of coupling for most protein 
machines is lower or equal to unity. Simultaneously, 
most protein enzymes display the Michaelis-Menten de- 
pendence of the asymptotic fluxes on the substrate con- 
centration. Gating the reactions by single pairs of con- 
formational transition substates is a necessary condi- 
tion for the conformationally fluctuating enzymes to 
obey the Michaelis-Menten kinetics [1, Hg]. There are 
thus solid grounds to suppose that the theory presented 
in Section III is applicable in the description of ac- 
tion of most biological machines. Doubts can be set- 
tled by analysis of time correlation functions of the di- 
chotomic noise observed in appropriate single-molecule 
experiments ^55^, '62j] . 

Of course, networks with gates comprising single tran- 
sition substates should be treated only as effective ones. 
The actual networks of conformational transitions are 
certainly much more complex. Various networks and 
systems of gates lead to the same or similar values of 
the quantities tm(^o {I J'}) in the expressions for the 
slippage functions Wi{Ai). Similarly, various networks 
make identical predictions of the statistical properties of 
the dichotomous noise observed [13, |6l|. It is a task 
for theoreticians to propose an algorithm of constructing 
the minimum effective networks that interpret the flux- 
force characteristics of the particular classes of protein 
machines. 

We tried to justify a hypothesis that the protein con- 
formational transition networks, like higher level biolog- 
ical networks: the proteome and the metabolome, have 
evolved in a process of self-organized criticality. A pro- 
posal follows from it to adopt evolving scale-free trees 
for the universal models of the conformational transition 
networks in the biomolecular machines. We assumed that 
the free energy-donating reaction (usually, the ATP hy- 
drolysis) is singly gated and proceeds through the main 
hub. In fact, the dependence of both the input and the 
output fluxes on the ATP concentration found in our 
simulations is of the Michaelis-Menten form, what agrees 
with many experiments. The universality of the ATP hy- 
drolysis is to be confronted with the fact that the main 
hub is very stable and evolves slowly. On the other hand, 
nodes with low connectivity evolve faster and can be fit- 
ted evolutionary, being good candidates for, if need be, 
either single or multiple entrance gate of the free energy- 
accepting process. 
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